clear;
clc;
addpath(genpath('./util'))
%大尺度衰弱参数
D0=1;
C0=10^(-30/10);
exponent_BS_LIS=2.2;
exponent_LIS_USER=2.8;

%BS位置
x_BS=0;
y_BS=0;
z_BS=25;

%LIS位置
x_IRS=0;
y_IRS=10;
z_IRS=40;

%用户位置范围
x_left=-50;
x_right=50;
y_left=10;
y_right=110;

%噪声功率
sigma=10^(-80/10)*10^(-3);

%发送功率
power=10^(45/10)*10^(-3);

%传输信噪比
rho = power/sigma/10^(14);

%BS-LIS距离 为防止运算溢出，故乘以10^(7)
d_BS_LIS=sqrt((x_BS-x_IRS)^2+(y_BS-y_IRS)^2);
[pathGain_BS_LIS] = pathLOS(d_BS_LIS,D0,exponent_BS_LIS,C0)*10^(7);

K = 4;
M=8;
NUM_timesamples = 20;

e=10^(-3);
for ueIdx = 1:K
    tau(ueIdx)= 0; %rand;sqrt(0.01)
end



L = 10:20:90;         % SNR dB

almax = 1;
almin = 0;

alpha_up = almax; %/rho
alpha_low = almin;%/rho

d=10;
f=10;
c1=2;
c2=2;
NUM=70;
P=20;
result_PSO1=zeros(1,length(L));
result_MonteCarlo1=zeros(1,length(L));
result_PSO2=zeros(1,length(L));
result_MonteCarlo2=zeros(1,length(L));
result_PSO3=zeros(1,length(L));
result_MonteCarlo3=zeros(1,length(L));

for i=1:NUM_timesamples
    i
    
    %生成用户路径损失
    for k=1:K
        x_user=rand*(x_right-x_left)+x_left;
        y_user=rand*(y_right-y_left)+y_left;
        d_IRS_user=sqrt((x_user-x_IRS)^2+(y_user-y_IRS)^2);
        %为防止运算溢出，故乘以10^(7)
        pathGain(k,1)=pathLOS(d_IRS_user,D0,exponent_LIS_USER,C0)*10^(7);
    end
    
    for j = 1:length(L)
        l=L(j);
        
        H1=randn(l,M)+j*randn(l,M);
        Hbar2 =randn(K,l)+j*randn(K,l);
        
        %生成信道
       [H1_1,Hbar2_1,T2_1] = getChannel2(l,M,K,0,H1,Hbar2,pathGain,pathGain_BS_LIS);
       [H1_2,Hbar2_2,T2_2] = getChannel2(l,M,K,1,H1,Hbar2,pathGain,pathGain_BS_LIS);
       [H1_3,Hbar2_3,T2_3] = getChannel2(l,M,K,10,H1,Hbar2,pathGain,pathGain_BS_LIS);
       
       
       [ F_best_PSO1,max_PSO1,alpha1 ]=AO(P,c1,c2,l,NUM,K,T2_1,Hbar2_1,tau,rho,M,H1_1,e,"PSO",alpha_low,alpha_up);
       sum_mon1=Sumrate_MonteCarlo(K,l,M,rho,alpha1,1000,T2_1,Hbar2_1,F_best_PSO1,H1_1,tau);
       
       [ F_best_PSO2,max_PSO2,alpha2 ]=AO(P,c1,c2,l,NUM,K,T2_2,Hbar2_2,tau,rho,M,H1_2,e,"PSO",alpha_low,alpha_up);
       sum_mon2=Sumrate_MonteCarlo(K,l,M,rho,alpha2,1000,T2_2,Hbar2_2,F_best_PSO2,H1_2,tau);
       
       [ F_best_PSO3,max_PSO3,alpha3 ]=AO(P,c1,c2,l,NUM,K,T2_3,Hbar2_3,tau,rho,M,H1_3,e,"PSO",alpha_low,alpha_up);
       sum_mon3=Sumrate_MonteCarlo(K,l,M,rho,alpha3,1000,T2_3,Hbar2_3,F_best_PSO3,H1_3,tau);
       
       result_PSO1(j)=result_PSO1(j)+max_PSO1;
       result_MonteCarlo1(j)=result_MonteCarlo1(j)+sum_mon1;
       
       result_PSO2(j)=result_PSO2(j)+max_PSO2;
       result_MonteCarlo2(j)=result_MonteCarlo2(j)+sum_mon2;
       
       result_PSO3(j)=result_PSO3(j)+max_PSO3;
       result_MonteCarlo3(j)=result_MonteCarlo3(j)+sum_mon3;

    end
end
plot(L,result_PSO1./NUM_timesamples,'r-o','LineWidth',1.5,'MarkerSize',6);
hold on;
plot(L,result_MonteCarlo1./NUM_timesamples,'b-*','LineWidth',1.5,'MarkerSize',6);
hold on;
plot(L,result_PSO2./NUM_timesamples,'r-o','LineWidth',1.5,'MarkerSize',6);
hold on;
plot(L,result_MonteCarlo2./NUM_timesamples,'b-*','LineWidth',1.5,'MarkerSize',6);
hold on;
plot(L,result_PSO3./NUM_timesamples,'r-o','LineWidth',1.5,'MarkerSize',6);
hold on;
plot(L,result_MonteCarlo3./NUM_timesamples,'b-*','LineWidth',1.5,'MarkerSize',6);
grid on;

P1=result_PSO1./NUM_timesamples
M1=result_MonteCarlo1./NUM_timesamples
P2=result_PSO2./NUM_timesamples
M2=result_MonteCarlo2./NUM_timesamples
P3=result_PSO3./NUM_timesamples
M3=result_MonteCarlo3./NUM_timesamples

h1=xlabel('L');
h2=ylabel('ergodic sum rate (bps/Hz)');
legend('Large-System Approximation Sum Rate (bps/Hz)','Ergodic Sum Rate (bps/Hz)');
